% make_river_cilmatology.m  3/17/2014  Parker MacCready
%
% makes river climatology for the rivers in the rtools database

% ideally we should do this from the Ecology database

clear

fn = 'riverFlow_1998_2009.mat';

r = load(fn);

[nt,nriv] = size(r.Qr_flow);

Qclim = nan(366,nriv);
Tclim = nan(366,nriv);

for yd = 1:366
    mask = r.Qr_yearday == yd;
    this_Q = r.Qr_flow(mask,:);
    this_T = r.T_riv(mask,:);
    Qclim(yd,:) = nanmean(this_Q);
    Tclim(yd,:) = nanmean(this_T);
end

% fill in last day, for Leap Years
Qclim(366,:) = Qclim(365,:);
Tclim(366,:) = Tclim(365,:);

if ( sum(isnan(Qclim(:))) + sum(isnan(Tclim(:))) ) > 0
    disp('WARNING: NaNs in climatology!')
end

for ii = 1:length(r.rname)
    rn = lower(r.rname{ii});
    % fix names
    if strcmp(rn,'duwamish'); rn = 'green'; end;
    if strcmp(rn,'hammahamma'); rn = 'hamma'; end;
    riv_name{ii,1} = rn;
end

yearday = [1:366]';

save riverFlow_clim Qclim Tclim riv_name yearday

if 0
    figure
    %
    subplot(211)
    plot(yearday,Qclim)
    ylabel('Qr (m3 s-1)')
    xlim([0,366])
    %
    subplot(212)
    plot(yearday,Tclim)
    ylabel('T (degC)')
    xlabel('Yearday')
    xlim([0,366])
end

% save csv files to load in python (cludgey)
%
% start by writing out the data
csvwrite('Qclim.dat',[yearday, Qclim]);
% then make a header line
name_line = sprintf('%s,',riv_name{:});
% strip the last comma
name_line = name_line(1:end-1);
% add the index name
name_line = ['yearday,',name_line];
% write the results to a csv file
fid = fopen('Qclim.dat','r');
fid2 = fopen('Qclim.csv','w+');
fprintf(fid2,'%s\n',name_line);
while 1
    tline = fgetl(fid);
    if ~ischar(tline), break, end
    fprintf(fid2,'%s\n',tline);
end
fclose(fid);
fclose(fid2);

% same thing for the temperature climatology
%
% start by writing out the data
csvwrite('Tclim.dat',[yearday, Tclim]);
% then make a header line
name_line = sprintf('%s,',riv_name{:});
% strip the last comma
name_line = name_line(1:end-1);
% add the index name
name_line = ['yearday,',name_line];
% write the results to a csv file
fid = fopen('Tclim.dat','r');
fid2 = fopen('Tclim.csv','w+');
fprintf(fid2,'%s\n',name_line);
while 1
    tline = fgetl(fid);
    if ~ischar(tline), break, end
    fprintf(fid2,'%s\n',tline);
end
fclose(fid);
fclose(fid2);



